curev Subroutine

      pure subroutine curev(idim,t,n,c,nc,k,u,m,x,mx,ier)

      !  calling sequence:
      !     call curev(idim,t,n,c,nc,k,u,m,x,mx,ier)
      !
      !  input parameters:
      !    idim : integer, giving the dimension of the spline curve.
      !    t    : array,length n, which contains the position of the knots.
      !    n    : integer, giving the total number of knots of s(u).
      !    c    : array,length nc, which contains the b-spline coefficients.
      !    nc   : integer, giving the total number of coefficients of s(u).
      !    k    : integer, giving the degree of s(u).
      !    u    : array,length m, which contains the points where s(u) must be evaluated.
      !    m    : integer, giving the number of points where s(u) must be evaluated.
      !    mx   : integer, giving the dimension of the array x. mx >= m*idim
      !
      !  output parameters:
      !    x    : array,length mx,giving the value of s(u) at the different points. x(idim*(i-1)+j) will
      !           contain the j-th coordinate of the i-th point on the curve.
      !    ier  : error flag
      !      ier = 0 : normal return
      !      ier =10 : invalid input data (see restrictions)
      !
      !  restrictions:
      !    m >= 1
      !    mx >= m*idim
      !    t(k+1) <= u(i) <= u(i+1) <= t(n-k) , i=1,2,...,m-1.
      !
      !  other subroutines required: fpbspl.
      !
      !  references :
      !    de boor c : on calculating with b-splines, j. approximation theory 6 (1972) 50-62.
      !    cox m.g.  : the numerical evaluation of b-splines, j. inst. maths applics 10 (1972) 134-149.
      !    dierckx p. : curve and surface fitting with splines, monographs on numerical analysis, oxford
      !                 university press, 1993.
      !
      !  author :
      !    p.dierckx
      !    dept. computer science, k.u.leuven
      !    celestijnenlaan 200a, b-3001 heverlee, belgium.
      !    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
      !
      !  ..scalar arguments..
      integer, intent(in) :: idim,n,nc,k,m,mx
      integer, intent(out) :: ier
      !  ..array arguments..
      real(RKIND), intent(in) :: t(n),c(nc),u(m)
      real(RKIND), intent(out) :: x(idim,m) ! x has size (mx), assume 2d (idim,m)
      !  ..local scalars..
      integer :: i,j1,k1,l,ll,l1,nk1
      real(RKIND) :: arg,tb,te
      !  ..local array..
      real(RKIND) h(MAX_ORDER+1)
      !  ..
      !  before starting computations a data check is made. if the input data
      !  are invalid control is immediately repassed to the calling program.
      ier = FITPACK_INPUT_ERROR
      if (m<1) return

      ! Check monotonic
      if (m>1 .and. any(u(2:m)<u(1:m-1))) return

      ! Check enough output space
      if (mx<(m*idim)) return

      ier = FITPACK_OK

      !  fetch tb and te, the boundaries of the approximation interval.
      k1  = k+1
      nk1 = n-k1
      tb  = t(k1)
      te  = t(nk1+1)
      l   = k1
      l1  = l+1
      !  main loop for the different points.
      eval_points: do i=1,m

        ! fetch a new u-value arg.
        arg = min(max(u(i),tb),te)

        ! search for knot interval t(l) <= arg < t(l+1)
        do while (.not.(arg<t(l1) .or. l==nk1))
          l = l1
          l1 = l+1
        end do

        ! evaluate the non-zero b-splines at arg.
        h = fpbspl(t,n,k,arg,l)

        ! find the value of s(u) at u=arg.
        ll = l-k1
        do j1=1,idim
          x(j1,i) = dot_product(h(1:k1),c(ll+1:ll+k1))
          ll = ll+n
        end do
      end do eval_points

      return
      end subroutine curev